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ABSTRACT 

Radiative pressure exerted by line interactions is a prominent driver of outflows in astrophysi- 
cal systems, being at work in the outflows emerging from hot stars or from the accretion discs 
of cataclysmic variables, massive young stars and active galactic nuclei. In this work, a new 
radiation hydrodynamical approach to model line-driven hot-star winds is presented. By cou¬ 
pling a Monte Carlo radiative transfer scheme with a finite-volume fluid dynamical method, 
line-driven mass outflows may be modelled self-consistently, benefiting from the advantages 
of Monte Carlo techniques in treating multi-line effects, such as multiple scatterings, and 
in dealing with arbitrary multidimensional configurations. In this work, we introduce our ap¬ 
proach in detail by highlighting the key numerical techniques and verifying their operation in a 
number of simplified applications, specifically in a series of self-consistent, one-dimensional, 
Sobolev-type, hot-star wind calculations. The utility and accuracy of our approach is demon¬ 
strated by comparing the obtained results with the predictions of various formulations of the 
so-called Cak theory and by confronting the calculations with modern sophisticated tech¬ 
niques of predicting the wind structure. Using these calculations, we also point out some 
useful diagnostic capabilities our approach provides. Finally we discuss some of the current 
limitations of our method, some possible extensions and potential future applications. 
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1 INTRODUCTION 

Predicting massive star evolution is dramatically complicated by 
the presence of powerful winds, constituting an importan t mass- 
loss mechanism (see, for example, overview in iKudritzki & Pulsl 
2000). The evolutionary timescales of the star and its luminosity 
are, for example, significantly affected by this continuous mass 
loss. In addition to influencing the stars themselves, these winds 
also affect their environments by injecting energy and momentum 
into the interstellar medium and supplying chemically enriched ma¬ 
terial. For a detailed study of all these processes and effects a firm 
understanding of the mechanisms driving stellar winds is required. 

In the case of hot O and B stars, the main driver for the mass 
outflow has been identified as the momentum transfer mediated 
by a multitude of at omic line interactions of the radiation fiel d 
in the wind material l lLucv & Solomonlll97d : ICastor et alJI 1975h . 
Over the years, this picture has been significan tly refin ed and a 
“standard model” (see, for example, review by IPuls et al.11200a) 
for li ne-driven winds ha s emerged, with i mport a nt contributions 
from Castor et~ail il975t) : IPauldrach et alJ dl986h : IKudritzki et alJ 
(1989). With the inclusion of multi-line effects, in particular by 
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IVink et ah , ( 2000 ). this model has become very successful in de¬ 
scribing hot-star winds and predicting their properties. However, 
observational diagnostics have identified a number of puzzles 
that challenge this “standard model” (see again IPuls et al.ll2008l) . 
Among these, the so-called clumping problem takes a prominent 
role. Contrary to the assumptions of the “standard model”, hot- 
star winds seem to exhibit strong temporal and spatial variabil¬ 
ity. Intense effort has been invested into understanding the prop¬ 
erties and the nature of these clumpy outflows, both from an ob¬ 
servational and theoretical viewpoint. In light of these efforts, we 
have developed a new radiation hydrodynamical approach to nu¬ 
merically model line-driven mass outflows self-consistently. Due 
to its reliance on Monte Carlo radiative transfer techniques, this ap¬ 
proach should be well-suited for time-dependent and multidimen¬ 
sional self-consistent studies of hot-star winds. In previous works, 
a Monte Carlo radiative transfer scheme has already been suc¬ 
cessfully coupled with a fluid dynamical method, resulting in the 
radiation hydro dynamical approach Mcrh jNoebauer et alJl2012l : 
Noebauer 2014). In this study we present extensions to MCRH and 
demonstrate the utility of this method for the investigation of line- 
driven winds. Thi_s_work has been partially carried out during the 
PhD project of lNoebauetl J20l4 ). 

The presentation of our approach is structured as follows: in 
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2 Noebauer et al. 


Section [2] we briefly recapitulate some basic concepts of hot-star 
winds and outline different solution strategies to address this prob¬ 
lem. We also highlight the advantages of Monte Carlo radiative 
transfer techniques for the study of line-driven mass outflows. We 
continue with a detailed description of the numerical techniques 
relevant for our approach in Section [3] In Section [4] we present a 
series of simple wind simulations to test our approach and in Sec¬ 
tion [5] we confront our method with an alternative, modern and 
sophisticated technique of predicting the wind structure. The ob¬ 
tained results and potential future improvements of our method are 
discussed in detail in Section[6] 


2 T HE ORY 

2.1 Solving the stationary hot-star wind problem 


Numerous techniques and approaches have been used to solve for 
the structure of hot-star line-driven winds and derive their prin¬ 
ciple properties. Many of these share the common feature of ad¬ 
dressing the problem in a “radiation-hydrostatics” fashion. By as¬ 
suming that the mass outflow is smooth, spherically-symmetric 
and in a stationary state, the wind structure may be deriv ed from 
globa l energy conservation considerati ons (e.g. lAbbott & Lucvl 
1 19851 iLucv & Abbottl fl993l. H l2004h o r by solving the mo- 
mentum equation alone (e. g. Castor et~ail 1 19751 . IPauldrach et alJ 
1 1986L [Kudritzki et al J1 1 98 c l) . If included, inhomogeneities in the 
mass outflow, i.e. clumps, are typically incorporated in a parame- 
terised fashion , for example, by introducing a clumping factor (e.g. 
[Schmutzil 1997l ; fl7amann & Koesterkeil 19981) . In the simplified, sta¬ 
tionary situation, the momentum equation reduces to 


dtt 1 d P GM * 

U ~A I A ^ 2 — 9rad, 

dr p dr r* 


( 1 ) 


after introducing the wind velocity u, its density p, the thermo¬ 
dynamic pressure P, the stellar mass M* and the radiative accel¬ 
eration p ra d. Solving this equation is still very challenging, since 
the bulk of the radiative acceleration derives from the momentum 
transfer mediated by interactions of the radiation field emitted by 
the central star with a large number of atomic line transitions. The 
contribution due to Thomson scattering may be taken into account 
by reducing the stellar mass in the gravity term by a factor (1 — T e ), 
using the Eddington factor with respect to electron scattering, T e . 

Various approaches have been used to calculate the line¬ 
driving force. We briefly sketch the basic prin ciple of t he so -called 
Castor-Ab bott-Klein (Cak) theory (after Castor et alJ 197.4. see in 
particular lAbbottl 1982 . Pauldrach et al.l 1986L Kudritzki et all 19891 
for improvements of the original approach) since it will be used for 
comparison in this work. lCastor et alJdl975h found an approximate 
analytic expression for the line-driving force. In particular, the line¬ 
driving force is expressed as a multiple of the electron-scattering 
contribution and this so-called force-multiplier is approximated by 
a power-law parameterisation. Given this form of the radiative ac¬ 
celeration, an analytic solution to the momentum equation may be 
identified. By performing a critical point analysis, the wind velocity 
is found to follow a /3-type law 


/ R 

u(r) = Uoo f 1- ^-J (2) 

with /3 = 0.5 in the original CAK case of ICastor et al .1 dl975t) . in 
which the central star is assumed to radiate like a point source and 
with /3 ~ 0.8 if the finite extent of the star is t aken into account 
(see lPauldrach et alii 19861 : iKudritzki et alJI 19891) . Throughout this 


work, we reference this improved formulation of the CAK approach 
as Mcak (“Modified-CAK”). The mass-loss rate of the wind scales 
as 

M oc C(a)(kL*)“ [M*(l — r e )] 1_ “ . (3) 


In addition to the stellar parameters (luminosity L*, photospheric 
radius 7?* and central mass M*), the parameters of the power-law 
fit, k and a. appear. The functional form of the coefficient C(a) dif¬ 
fers slightly between the CAK and Mcak approaches. More details 
about C(a) and the Cak/Mcak methods, as used in this work, are 
provided in Appendix lAl 

A significant advantage of these radiation-hydrostatic tech¬ 
niques lies in the analytic expressions for the wind structure they 
either provide or in the po ssibility to efficiently solve the prob¬ 
lem numerically (see, e.g., Abbott & Lucv 1985:j; IPauldrach et alJ 
1 1994l:lvink et al.ll999l : IPuls et alj2005l) . even if more sophisticated 
treatments of the line-driving force are included. In both cases, 
wind models may be constructed for entire grids of different stel¬ 
lar parameters (e.g.|Abbottll982l : IVink et alj200d : IPuls et alj2005l : 
12 ). 


IMuiires et alJ 120121) . A drawback lies in the omission of a self- 

consistent treatment of temporal and spatial variations in the wind 
flow. 


2.2 A dynamical approach 

To account for, follow and solve the temporal and spatial vari¬ 
ability of stellar winds, the full radiation hydrodynamical problem 
the r adiatively-driven ma s s outflow constitutes has to be addressed 


(e.g. Owocki et al. 1988; Feldmeier 1995; Feldmeier et al 


llPessart & Qwockill2003l; D essar 120041; iDessart & Owockiil 


19971; 


20051) . 


In such an approach, a numerical solution to the full radiation hy¬ 
drodynamical equations describing the conservation of mass, mo¬ 
mentum and energy has to be formulated. In contrast to a pure 
hydrodynamical problem, the energy and momentum transfer be¬ 
tween the radiation field and the fluid material are included. These 
terms are found by solving the co-evolution of the radiation field 
in parallel. However, due to the complexity of the full radiation hy¬ 
drodynamical problem, a simplified treatment of the radiation field 
and radiative transfer is typically adopted (c.f. lOwocki et alJI 19881 : 
iFeldmeierll 1 995l ; |Pessart & Owockill2003l) . The radiation hydrody¬ 
namical equations are given in Appendix Q3] in the form in which 
they are used in this study. 

In this work we propose a new numerical approach to per¬ 
form such dynamical self-consistent simulations of hot-star winds. 
The distinguishing feature of our method is the incorporation of 
Monte Carlo radiative transfer techniques. In particular, we build 
upon the general-purpose Mont e Carlo radiati o n hyd r odynamics 
schem e MCRH, developed by iNoebauer et alJ d2012h : iNoebaueil 
d2014l) and adapt it to the line-driven wind environment. Other 
Monte Carlo-base d radiation hyd r odynamical techniques hav e 
been developed by Na yakshin et al.l ((2009), Acre man et al.l d2010h. 
lHaworth & Harries 1 20121) and most recently by Roth & Kasenl 
d2015l) . who place a particular e mphasis on the incorporati on of 
implicit Monte Carlo methods (see lFleck & Cummingslll97lh . and 
bv lHarriesI d2015h . 


2.3 Advantages of Monte Carlo techniques 

Monte Carlo techniques have become established as a competitive 
and successful numerical approach to solve radiative transfer prob¬ 
lems due to the increasing availability of computational resources 
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Line-driven winds with Monte Carlo radiation hydrodynamics 3 


and the more and more levels of sophistication added to the for¬ 
malism (see specifically, Abbott_&Lucy j_985L iLucvll 1999al iLucvl 


lev I TO. 

I2002L lLucvll2003l . ICarciofi & Biorkmanl 20061) . The appeal of us 
ing Monte Carlo techniques derives from t heir advantages in treat 


ing interaction physics (see, for example, IPozdnvak ov et alj ^983 
for M ont e Carlo calculations of comptonisation or iKasen et al 
120061 and lKromer & Simll2009l . for presentations of fully-fledged 
three-dimensional Monte Carlo codes to solve multi-frequency ra¬ 
diative transfer in Type la supernova ejecta) and from the ease 
with which arbitrary geom etrical configurations are treated (see, 
for example, Camps et al.ll2013L describing a Monte Carlo radia¬ 
tive transfer scheme on a Voronoi mesh). These advantages all de¬ 
rive from the microscopic viewpoint that Monte Carlo techniques 
take by solving the propagation history of a large number of rep¬ 
resentative photons probabilistically. They are very relevant for the 
study of hot-star winds, in which a multitude of atomic line inter¬ 
actions are the main driver, multiple scatterings frequently occur 
and a clumpy irregular structure of the mass outflow is expected. 
Thus, Monte Carlo techniques have been frequently used to ad- 
dr ess the radiative transfer problem in hot - star w i nds, for example 
bv I Abbott & Lucvl Jl985l).lL ucv & A bbotdi 1 993ll. lSchmutzl J 1 9971) 


IVink et al] dl999h. Ivink et al] OOOOh. I Si nil J 20G4|) 7 Tm u 11 er & Vinkl 

1 2008h . Muiires et al. 1 2012h . lSurlan et al d2012h . 

In this work, we aim at exploiting the advantages of Monte 
Carlo techniques for solving the line-driving problem in a dynam¬ 
ical and self-consistent approach by using a Monte Carlo radiation 
hydrodynamical technique. This strategy distinguishes itself form 
previous Monte Carlo-based investigations by its treatment of the 
coupled co-evolution of the radiation-wind state. Moreover, the ra¬ 
diative acceleration is determined locally, which sets us apart from 
some e arlier Monte Carlo-ba s ed studies, such a s |Abbott&Lucy| 
dl985h . lLucv & Abbott! Jl993h . lVink et all (1 19991) and lsinj i |2004) . 


in which the momentum transfer onto the wind material was de¬ 
rived from global considerations without insisting on local con¬ 
sistency (however, see lMiiller & Vinkll2008l . for developments to¬ 
wards a consistent local force balance in Monte Carlo-based ap¬ 
proaches). 


3 NUMERICAL METHODS 

In order to adequately address the radiation-matter coupling prob¬ 
lem in the stellar wind environment, we have extended our numeri¬ 
cal framework whose underlying m ethodolog y and first implemen¬ 
tation has been described in lNoebauer et al.Hl2012l) . In the follow¬ 
ing, the basic principles of our approach are briefly recalled and the 
new extensions described in detail. 


3.1 Basic operator-split approach 

In our approach, a simple operator-splitting scheme is employed, 
coupling a time-dependent Monte Carlo radiative transfer simula¬ 
tion with a finite-volume fluid dynami cal calculation. In the latter , 
the piecewise parabolic method (PPM. IColella & Woodwardl 19841) 
is used to construct a series of Riemann problems at the interfaces 
of the computational do main, which are solved with a standa rd Rie¬ 
mann solver (HLLC, see lToro et al J 1 994l : lBatten et alll997l) to de¬ 
termine the mass, energy and momentum flux through the inter¬ 
faces and eventually the evolution of these conserved quantities in 
the cells. In the following radiative transfer step, a Monte Carlo 
simulation is performed to determine the evolution of the radiation 
field and reconstruct the energy and momentum transfer between 


the radiation field and the fluid material. In the final step of the 
splitting approach, the momentum and energy of the fluid are al¬ 
tered in accordance with these transfer terms. 

This basic outline of the numerical scheme remains valid for 
the stellar wind applications. Alterations and extensions of the 
scheme only concern the inclusion of additional physical effects or 
adopted simplifications that are relevant for an adequate treatment 
of the stellar wind environment. In particular, the gravitational field 
originating from the central star is taken into account, an isother¬ 
mal treatment of fluid dynamics is adopted and the radiative trans¬ 
fer scheme is generalised to incorporate frequency-dependent res¬ 
onant line interactions. These alterations are systematically intro¬ 
duced and presented in detail in the following sections. 


3.2 Isothermal hydrodynamics and gravity 


In this work, we reduce the complexity of the hydrodynamical cal¬ 
culations in our modelling approach of hot-star winds by adopting 
an isothermal approximation. This simplification, which has been 
used in numerous investigations of line-driven hot-star winds (e.g. 


1:1 Abbott & Lucvl[l985l: lOwocki et alJI 1988c Ivink et aD 


199 fe00~diSinj2004l ). is justified by the continuous irradiation of 

the wind material by the central star, which together with a char¬ 
acteristic cooling time scale that is faster than the flow time scale 
sustains the wind material roughly at the effective temperature of 
the star dKlein & CastoiilT978b . In the isothermal approximation, 
only the mass and momentum conservation equations have to be 
addressed in the fluid dynamical calculation. With the equation of 
state of an isothermal flow directly connecting the fluid density and 
pressure 


T) 2 


(4) 


a solution of the energy conservation equation is not required. In 
the case of an ideal gas with constant temperature Tl so and mean 
molecular weight g,, the isothermal sound speed ai so is 


CLiso — 


/ 


ksTiso 

P 


(5) 


We modify the fluid dynamical solution procedure in the 
isothermal version of Mcrh by reconstructing only the primitive 
variables density and the velocity. Using the isothermal equation 
of state, the pressure may be directly determined and the Riemann 
problems at the cell interfaces solved. Here, we use a particular 
version of the HLLC Riemann solver, adopted from the Athena 
co dJ3 dStone et alil2008 ). which is tailored to the particular form 
of the simple waves in isothermal flows. 

An important contribution to the momentum balance in hot- 
star winds derives from the gravitational pull exerted by the central 
object 


9 = 



( 6 ) 


Following IColella & Woodward! d 19841 ). we incorporate the effect 
of the star’s gravitational pull into the isothermal fluid dynami¬ 
cal calculation by accounting for the additional momentum density 
contribution 


A n ^ A n n . n+1 n + l\ , , 

Ap = -At(p g +p g ), (7) 


1 The source code of this finite-volume magneto-hydrodynamical code 
may be obtained from https : //trac. princeton. edu/Athena/ 
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4 Noebauer et al. 


in the final determination of the new fluid dynamical state at the 
end of the splitting step. Here, the superscripts n and n + 1 mark 
quantities at the beginning and the end of a time step with length 

At = t n+1 - r. 


3.3 Monte Carlo principles 


take convenient forms. Transformation rules relate quantities be¬ 
tween the se two frames^ These are derived in detail, for ex¬ 
ample, in iMihalas & Mihalasl d 1984tf and are listed partially in 
iNoebauer et alj 1 20121) . 

During their propagation, packets interact whenever they have 
covered a pathlength equivalent to the optical depth 


t = — In z, 


( 11 ) 


At the heart of the Monte Carlo radiative transfer machinery lies the 
discretization of the radiation field into a large number of quanta, 
so called packets, which are launched into the computational do¬ 
main. In a series of random number experiments, the packets de¬ 
scribe radiation-matter interactions and the temporal evolution of 
the radiation field. The basic layo ut of our Monte Carlo method has 
already been presented in detail in INoebauer et alj J20 1 2l ). Here, we 
only recall some important aspects and dedicate the bulk of the pre¬ 
sentation to the additions related to the frequency-dependent calcu¬ 
lation. 


3.4 Discretization, initialization and propagation 

We ad opt t he indivisible energy packet scheme of I Abbott & Lucvl 
(1985) and iLucvl l l 1 999,1) . Thus, Monte Carlo packets are inter¬ 
preted primarily as parcels of radiative energy. To model the con¬ 
tinuous illumination of the outflowing material by the central star, 
a number of packets are injected into the computational domain 
through the lower boundary in each radiative transfer step. Assum¬ 
ing that the star radiates as a black body with effective temperature 
T e s, 


N, 


inject 


477 Rl (7 R T e 4 ff At 
£ 


(8) 


new packets, each with energy e, are launched from the inner 
boundary (located at .R*) during a time step of length At. Neglect¬ 
ing any limb-darkening effect, the initial propagation direction of 
these packets follows 


P=Vz. 


(9) 


Here and in the following, p denotes the cosine of the directional 
angle between the trajectory and the symmetry axis and s repre¬ 
sents a random number, uniformly drawn from the interval [0,1], 
The frequency assignment of the packets is chosen to reflect the 
Planck functior0 


B(u,T) = 


2 hu 3 1 

c 2 ghu/k B T _ ^ ' 


( 10 ) 


This frequency assignment is adopted for this current study for sim¬ 
plicity and convenience. Sampling a realistic atmosphere model, 
however, does n ot pose a con ceptual or implementation challenge 
(see, e.g., I Abbott & Lucvin~985l) . 

After determining the initial properties of the packets, they 
are Jaunched and initiate the propagation process. As detailed in 
INoebauer et aD we follow a mixed-frame approach and 

perform the packet propagation procedure in the lab frame ( LF ) 
in which the computational grid is defined. However, all interac¬ 
tion processes are treated in the local co-moving frame ( CMF ), 
in which the fluid is at rest and where the material functions 


2 See for example Carter & Cashwell (1975), for an efficient algorithm to 
sample this function. 


which is determined for every packet in a random number experi¬ 
ment. The procedure with which the nature of the interaction event 
is determined is presented in the next section. Once the details of 
the interaction process are established, the packet properties are up¬ 
dated accordingly. We highlight the case of line interactions, which 
we treat as reson ant scatterings in the Sobolev approximation (af- 
ter lSob olev 1960; see also, for example. lLamers & Cassinel fill 19991 
for a summary of this approximation and Pauldr ach et alJI 19 86 for 
a discussion of its applicability to hot-star winds). In this case, the 
packet is assigned a new propagation direction in the CMF accord¬ 
ing to the Sobolev escape probability 


1 — exp(—T a ) 

p(m) oc -. (12) 

T h 

Here, r s denotes the Sobolev optical depth, which will be explicitly 
introduced in Section lXTl and which depends on the direction of the 
photon trajectory. The changes in the LF frequency and energy of 
the packet follow from the Doppler effect and from energy con¬ 
servation in the CMF. Assuming that the process of sampling the 
Sobolev escape probability returned the CMF propagation direc¬ 
tion po and accounting for the appropriate transformation laws, the 
resonant line scattering process results in the new packet quantities: 


a _ Mo + P 

M “ 1 + £m§ ’ 

e a = 7 2 (1 - /3p h )(l + /3po)e h , 
v 3 - = 7 2 (l - /3/z b )(l + /3po)v h . 


(13) 

(14) 

(15) 


Here and in the following, a subscribed “0” denotes CMF quanti¬ 
ties and the superscripts b and a distinguish packet properties be¬ 
fore and after the interaction process. Moreover, the standard no¬ 
tation of special relativity is adopted, introducing (3 = u/c and 
7 = 1/ -y/l — /3 2 . Above and in the following we keep full account 
of the relativistic terms and only in the final expression reduce the 
accuracy to 0(u/c), which is appropriate for modelling hot-star 
winds. 

In addition to the physical radiation-matter interactions, pack¬ 
ets also experience so-called numerical events due to the spatial 
and temporal discretization of the computational domain. When a 
packet intercepts a grid cell boundary, some packet properties, such 
as the optical depth distance to the next interaction, have to be re¬ 
determined to be compatible with the physical conditions in the 
target cell. At the end of the time step, the trajectories of all re¬ 
maining packets are interrupted and their properties stored to allow 
them to resume their propagation during the next simulation cycle. 


3.5 Optical depth summation 

The sole purpose of the packet propagation process lies in the de¬ 
termination of packet interaction histories from which the radiation 
field characteristics can be reconstructed. For the stellar wind ap¬ 
plication, this concerns primarily the line-driving force. Crucial for 
the determination of the packet trajectories is the decision about 
where and how the packets interact. As already mentioned, pack¬ 
ets experience physical radiation-matter interactions after having 
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accumulated an integrated optical depth equal to the random num¬ 
ber experiment outcome CD}- For the current work, the wealth of 
possible interaction processes has been restricted to include only 
frequency-dependent bound-bound processes, which we treat as 
resonant scatterings in the Sobolev limit (see previous section). All 
continuum processes are neglected, apart from Thomson scatter¬ 
ing, which is incorporated approximately by reducing the mass of 
the central star as outlined in Section 12711 

Once a packet interacts, the location of this event in optical 
depth space has to be translated into a physical position within the 
computational domain. This procedure is trivial in the grey case, 
since optical depth and pathlength only differ by the opacity, a 
constant multiplicative factor in the CMF, but complicated in the 
presence of frequency-dependent proces ses. We adopt a simplified 
version of the technique of iMazzali & Lucvl d 1993h to locate the 
line-interaction events packets perform. On its trajectory, a packet 
propagates freely to the Sobolev point of the next line with which it 
comes into resonance. Each time such a resonance point is reached, 
the optical depth is incremented by the full line optical depth of 
the corresponding transition. The packet undergoes an interaction 
once the value drawn in G3 is surpassed by the optical depth ac¬ 
cumulated. If this occurs during the instantaneous increases at one 
of the Sobolev points, the packet undergoes a resonant line inter¬ 
action, otherwise it may leave the current grid cell uninterrupted. 
This procedure may be easily extended to include additional inter¬ 
action types, in particular frequency-independent processes, such 
as Thomson scatterings fsee lMazzali & Lucvll993l) . but for the cur¬ 
rent work we have omitted to do so. 


3.6 Monte Carlo estimators 

In reconstructing the radiation field characteristics from the ensem¬ 
ble of packet interaction histories, we f ollow the volume-averaged 
estimator approach proposed bv lLucvl dl999al) and refined bv lLucvl 
J2003L120051) . This formalism aims at reducing the statistical fluc¬ 
tuations inherent to the Monte Carlo approach by increasing the 
number of contributions to the packet census. For the case of 
frequency-independent processes being the only interaction chan- 
ne l, adequate estimator s have already bee n derived and present ed 
bv lNoebauer et alJ d2012l) and similarly bv lRoth & Kasenl J20I5I ). 

To determine the radiative acceleration due to spectral line in¬ 
teractions, we consider the momentum transferred in such an event. 
Assuming that these interactions occur as resonant scatterings, a 
packet transfers 

Np = £ — [/x b -7 2 (Mo + /9)(! - /V’)] (16) 

momentum onto the material. Estimators for the radiation force can 
be obtained by summing over the transfer terms of all interacting 
packets. To reduce the statistical fluctuations in these estimators 
we follow the suggestion of lLucvl d 1999bt) and include all packets 
that come into resonance with a line and weight their contributions 
with the corresponding interaction probability given by (1 — e~ Ts ). 
Taking the forward-backward symmetry of the re-emission into ac¬ 
count, thus cancelling all terms that are of odd power in fj,Q, the 
following estimator for the radiation force [c.f. Equation dB5H due 
to line interactions are obtained: 

GLe = AV^Af B 1 - e ~ Ts ^ ~ ?)■ ( 17 > 

Here, the volume of a grid cell, AV appears. The superscript b 
has been dropped and only terms of 0(u/c) have been retained. 


Using this estimator, the radiation force is calculated and the fluid 
momentum updated in the final splitting step. 


3.7 Ionization and level population 


The strength of the different line transitions, as encoded in the 
corresponding values for the Sobolev optical depth, t s , depends 
on the excitation and ionisation balance in the wind. Thus, the 
Monte Carlo radiative transfer routine requires a separate calcula¬ 
tion step which determines the ionisation and the level population 
in e ach cell of the compu t ational grid . In thi s fir st study, w e fol- 
low I Abbott & Lucvl d 1 9851) . [Vink et al J d 1 999t) and lSiml d2004l) and 
adopt an approximate non-LTE treatment. We stress, however, that 
nothing in our radiative transfer or hydrodynamics formalism re¬ 
quires these approximations. A full, non-LTE scheme for calculat¬ 
ing ionisation and excitatio n sta tes could be incorporated following 
the methods described bv lLuc 3 Hoi see Section 1631 . 

Following our simplified strategy, we determine the ionisation 
balance b y apply ing the modified nebular approximation (see, e.g. 
Mazzali & I .ucx l 199.il l. 

Nj + l = (NeY 1 
Nj \W) 


x [0 + W (1 - 0)] 



N j+1 N e 

Nj 


Tr 


(18) 


This expression relates the number densities of two successive ion 
stages, Nj, Aj+i with the electron number density N e . Compared 
to a pure LTE calculation based on the Saha equation, whose results 
are denoted by the asterisks, modifications due to the dilution of 
the radiation field and due to recombination effects are taken into 
account. As a consequence, the dilution factor W, the ratio of the 
electron and radiation temperatures, T e and Tr, and the fraction 
C,j of recombination processes going directly to the ground state 
appear in the above formulation. 

In the determination of the level population, we again follow 
lAbbott & Lucvl (Il985l) and identify levels with no permitted elec¬ 
tromagnetic dipole transitions to lower energy levels as metastable; 
these levels are assumed to obey Boltzmann statistics 



with “1” denoting the ground state. For all other levels, the effect 
of radiative processes on the level population is crudely taken into 
accou nt by including a dilution factor (see Ivink et alJ[l999l : ISiml 
1 20041) 



For the current work, we do not solve Equation d 1 8b iteratively 
to determine the ionisation state. Instead, we use a predefined char¬ 
acteristic electron dens ity N e /W and calculate the ionisation ac¬ 
cordingly (see, e.g. lAbbotdH982l . for a comparable strategy) in all 
our wind simulations. Also, the radiation temperature is treated in 
a simple fashion according to the prescription 


Tr = T e = Tefr 


( 21 ) 


rather than relying on predictions of realistic atmosphere models. 

With the excitation and ionisation balance accessible through 
Equations ( 118b . d 19b and ( 120b in each cell of the computational 
grid, the Sobolev optical depth of a line transition encountered on a 
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Figure 1. Illustration of the linear interpolation scheme to ensure continuity 
in the CMF packet frequency. The grid cell boundaries are treated as loca¬ 
tions where the fluid velocity and thus the CMF frequency is known at all 
times. The fluid flow speed there is obtained by linear interpolation (dashed 
line) between the volume-averaged values of the finite-volume scheme. The 
velocity in the grid cells is then obtained by another linear interpolation step 
between the boundary values (solid blue line). 

packet trajectory along the direction p may be calculated according 
to 


c 

(niP)r a 

(22) 

r “ Vo [(1 - 

m 2 )t + m 2 £L/ 

7re 2 

Kip = - JlTil 

(i _ Ih.EL] 

(23) 

m e c 

V ni g u ) 



Here, uq denotes the rest frame frequency of the transition and /; 
its oscillator strength. Also, the statistical weights gi iU associated 
with the lower and upper energy levels connected by the transition 
and the electron charge and mass, e and m e , appear. Finally, we 
stress again that the above expressions only depend on the phys¬ 
ical conditions at the Sobolev point, an essential feature which is 
highlighted by the subscribed r s . 

3.8 Velocity interpolation 

As photons propagate through the wind material, their CMF fre¬ 
quency is continuously Doppler-shifted by the varying wind veloc¬ 
ity. To reproduce this behaviour in the Monte Carlo simulation, an 
interpolation scheme is required to reconstruct the evolution of the 
wind velocity within the computational grid cells. For the current 
work, we do not expect the formation of shocks but the presence of 
a smoothly varying velocity field. Thus, in constructing the algo¬ 
rithm, an artificial introduction of discontinuous jumps in the inter¬ 
polated wind velocity should be avoided: otherwise, packets may 
skip frequency regions which may be populated by line transitions. 
We prevent this by devising an interpolation scheme that insists on 
continuity in the reconstructed wind velocity at the interfaces be¬ 
tween grid cells. First, the velocity at the interfaces are determined 
by linear interpolation between the cell-averaged values provided 
by the finite-volume fluid dynamical calculation. The interface val¬ 
ues are then used in a second linear interpolation step to determine 
the evolution of the wind velocity within the cell. Figure Q] graphi¬ 
cally illustrates this procedure. 

Similarly to the spatial considerations, care has be taken when 
packets reach the end of the time step. In order to ensure consis¬ 
tency in the CMF frequency between successive time steps, the ve¬ 
locity is also interpolated linearly in time between the results of two 


preceding fluid dynamical splitting steps. Otherwise the changes of 
the fluid velocity due to fluid dynamics and the radiation-matter 
coupling could introduce discontinuous jumps in the CMF fre¬ 
quency. 


3.9 Boundary conditions 

In the hot-star wind calculations, the inner boundary of the compu¬ 
tational domain is set to the stellar surface, f?*. So-called ghost 
cells, which are required for the fluid dynamical reconstruction 
step at the domain boundaries, are inserted below this radius. The 
fluid state in these inner ghost cells is found by using a tempo¬ 
rally constant hydrostatic density stratification and by extrapolat¬ 
ing the velocity field from the innermost region of the domain into 
the ghost cells using a low-order polynomial. For most wind cal¬ 
culations presented in this work, the hydrostatic profile has been 
normalised (arbitrarily) to po = 10 -11 gem -3 . However, no sig¬ 
nificant changes in the calculated wind structure was found if this 
value was modestly varied. In the future, predictions from atmo¬ 
spheric models may be used for this process instead. With the wind 
density fixed but the possibility of the velocity in the ghost cells 
to float, the mass flow through the boundary may quickly adjust it¬ 
self to the conditions in the wind. Thi s situ ation closely resembles 
the boundary conditions presented bv lOwocki et al.l ( il988l) for their 
time-dependent studies. 

As detailed in Section [3~4l the inner edge of the domain con¬ 
stitutes an inflow boundary for the radiation field. Any Monte Carlo 
packet that is back-scattered through the inner boundary is dis¬ 
carded from the ensuing propagation process. 

The outer edge of the computational domain is transparent to 
the radiation field. All escaping packets are recorded to construct 
spectra or other diagnostics during the postprocessing steps. Re¬ 
garding the fluid state, the outer boundary simulates the outflow of 
wind material. To this end, the values for the fluid variables in the 
ghost cells are found by extrapolating the wind flow. To minimise 
the effect of stochastic fluctuations introduced by the Monte Carlo 
radiative transfer step, a linear regression scheme is used. 

This outline of the boundary treatment concludes the descrip¬ 
tion of the numerical techniques which are relevant for this study 
and which have been added to Mcrh. The performance of all these 
extensions has been verified in a series of test calculations. A de¬ 
tailed description of this procedure and the different tes^problems 
used during this process may be found in lNoebauen ( 120141) . 


4 HOT-STAR WIND TEST CALCULATIONS 

We now aim at demonstrating the capability of our approach to 
self-consistently solve the line-driving problem in hot-star winds 
by considering an idealised representation of the problem. Specif¬ 
ically, we test the basic accuracy of our method by exploring 
whether it produces results compatible with the Cak and Mcak 
theory when comparable assumptions about the physical effects at 
work are adopted (smooth flow, Sobolev approximation, etc.). At 
the same time, we investigate how these simple calculations com¬ 
pare with simulations in which the full details as outlined in Section 
[3] are incorporated. Consequently, this series of wind simulations, 
constitutes the basic testing and validati on step before we confront 
our method with the technique of lMiiller & Vinkl d2008i) in Section 
[5] Finally, in Section[6] we discuss how some of the simplifications 
and approximations adopted in the current work may be relaxed and 
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Table 1. Overview of all chemical elements, that are taken into account in 
our wind calculations. For each species, the included ionisation stages are 
listed as well. 


Element 

Ions 

Element 

Ions 

H 

i,ii 

He 

mi 

C 

I IV 

N 

I-VI 

O 

I-VI 

F 

I-VI 

Ne 

I VI 

Na 

I-VI 

Mg 

I-VI 

A1 

I-VI 

Si 

I-VI 

P 

I-VI 

S 

I-VI 

Cl 

I-V 

Ar 

I-V 

K 

I-V 

Ca 

I-VI 

Ti 

I-VI 

Cr 

I-VI 

Mn 

I-VI 

Fe 

I-VI 

Co 

I-VI 

Ni 

I-VI 




Table 2. Properties of £-Puppis according to Puls et alj il996l) . We adopt 
these stellar parameters in our wind test calculations. 


Parameter 

Value 

Mi, 

52.5 M 0 

Li, 

10 6 Lq 

T e S 

4.2 X 10 4 K 


a more realistic description of the line-driving problem be achieved 
within our numerical framework in the future. 



Figure 2. Results of the direct determination of the force multiplier as a 
function of the optical depth parameter t, performed with MCRH (orange) 
according to Equation Cm). In the transition region, between the optically 
thick and thin regimes, a power-law fit according to Equation JaT} has been 
performed (dashed blue). The parameters of this fit are included in Table [3] 
and the predicted wind state (in terms of t ), according to the C AK theory and 
the obtained k and a values, is marked by a green cross. For comparison, the 
force multiplier values of lAbbottl ll982l) , determined for a star with T e ff — 
4 X 10 4 K and N e /W = 1.8 X 10 s , 1.8 X 10 11 , 1.8 X 10 14 cm“ 3 are 
included (grey open symbols). 


4.1 Atomic data 


4.3 Cak fitting 


In all wind calculations, we assume that the material has solar com¬ 
position. We include the elements listed in Table Q] and adopt the 
corresponding abundances of lAsplund et al.l 120091) . Table |T] also 
identifies the ionisation stages that are taken into account for each 
element. Information about the atomic structure and the line tran¬ 
sitions in these ions are taken from the sa me data base used in the 
study of colliding winds bv IParkin & Sim (2013), which is based 
on the Kurucz^ Bell| ( 119950 atomic data set. To reduce the com¬ 
putational effort in the Monte Carlo radiative transfer steps, we do 
not account for all line transitions recorded in the data set, but only 
those with log gf > —6. We have explicitly verified that the inclu¬ 
sion of more weak lines does not affect the outcome of our simula¬ 
tions. When required, the recombination fractions for the included 
ions, are adopted from the PYTHON code ( iLong & Kniggei 

l2002h . 

4.2 Stellar parameters 

We perform all wind calculations for fixed sets of stellar parame¬ 
ters. For the simulations series constituting the basic testing pro¬ 
cess, we consider a system that is similar to the well-studied O-star 
£-Puppis. The basic stellar parameters for these calculations are 
listed in Tableland are adopted from lPuls et al.l ( 119961) . These val¬ 
ues imply a stellar radius of f?* = 1.317 x 10 12 cm. Additionally, 
we use N e /W = 10 15 cm -3 in the £-Puppis calculations. Per¬ 
forming an ionisation and excitation calculation as outlined in the 
previous section, we find a mean electron scattering cross section 
of <7 e = 0.34cm 2 g _1 throughout the wind, which corresponds 
to the Eddington factor T e = 0.5. In all MCRH wind simulations 
presented in the following, we account for the effect of electron 
scattering by reducing the stellar mass by the factor (1 — T e ), as 
described in Sections [2. 1 l and [3~5l 


In the CAK theory, the wind structure may be readily determined 
from the stellar properties and the power-law parameters k and a 
of the force multiplier. However, the wind characteristics, mass- 
loss rate and terminal wind speed, are sensitive to the exact val¬ 
ues of these parameters. Given our aim at assessing the utility of 
our approach for solving the line-driving problem, we refrain from 
consul ting previous studies, such as lAbbottl dl982h or IPauldracil 
(1987), which provide these parameters for a wide range of stel¬ 
lar conditions. Instead, we use the Monte Carlo routine of MCRH 
to determine the values for k and a. This way, we ensure that 
the Mcrh-Cak/Mcak comparisons are not obscured by differ¬ 
ences in atomic data sets, the ionisation and excitation treatments 
or the stellar parameters. In particular, we determine the va lues for 
k and a by carrying out the force-multiplier summation (c.f. l Abbott] 
Il982h 


M(t) = £ 


F„ 0 Aurt 1 
F t 


1 — exp 


mt \ 

Az/ D CTe ef ) 


(24) 


explicitly with the MCRH modules for a large number of different 
values of the dimensionless optical depth and by performing the 
power-law fit according to Equation ( IA4b afterwards. The Doppler 
width, Aisd = voUt^/c, according to the thermal motion u t h [see 
Equation ( IA3H . is used together with a reference specific electron 
scattering cross section al et [see Equation ( fAP l. The results of the 
direct summation procedure are shown in Figure|2] together with the 
fitting curve. As a reference, the results of lAbbottl Jl982n . obtained 
for comparable physical conditions are included. In general, both 
calculations agree, but noticeable differences are present, which 
are most likely a consequence of our simplified treatment of ion¬ 
isation and excitation, differences in the atomic data sets and of the 
use of realistic stellar atmosphere models instead of a simple black 
body bv l Abbott! j 19821) . Given the purpose of our calculations, these 
differences are irrelevant but highlight the need for consistency be- 
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Table 3. Overview of the main Cak parameters determined in the fitting 
process and used throughout this work. In addition, the wind properties, as 
predicted by the Cak theory according to these parameters are listed in the 
bottom rows. The adopted values for N e /W, po and T e are included for 
completeness. 


CAK Parameter/Wind Property 

Value 

k 

0.381 

a 

0.595 

t 

0.035 

M(t) 

2.80 

M 

4.51 x 10 -5 Mq yr —1 

tioo 

881kms _1 

N e /W 

10 15 cm -3 

PO 

10~ n gem -3 

T e 

0.502 


tween the Mcrh and CAK calculations. The final CAK parameters, 
which are used throughout this work, are listed in Table[3] Based on 
the force multiplier parameters, the CAK predictions for the mass- 
loss rate and the terminal wind speed are calculated and also pro¬ 
vided in the table. Including the finite-cone effect according to the 
Mcak approach increases the terminal wind speed by a factor of 
~ 2.68 to about Uoo ~ 2360kms _1 and reduces the mass-loss 
rate to M ft! 2 x 10 -5 M 0 yr -1 . With these modifications, the ter¬ 
minal wind speed is compatible with previous investigations of the 
wind of ([-Puppis (e.g. lPuls et al.ll996f) . However, the density in our 
wind is too high. Specifically, the mass-loss rate is about a factor 
of 3.5 larger than established bv IPuls et al.l C 199fil ). Again, our sim¬ 
plified treatment of the wind ionisation and excitation conditions is 
most likely the cause for these discrepancies. We stress once more, 
however, that as long as we use the simplified description of the 
wind conditions consistently in the Cak/Mcak calculations and 
the Mcrh simulations, the conclusions drawn from comparing the 
corresponding results are unaffected. 


4.5 Simulation series layout 


Using the general parameters outlined in the previous sections, we 
perform a series of MCRH simulations of a (^-Puppis-like wind. In 
this series we successively increase the level of physical detail (see 
Section l4~7l l. For the first stage, the point-source approximation and 
an unattenuated radiation field will be used, similar to lCastor et alj 
Il975l) . In the second set of calculations, the finite extent of the 
star will be taken int o ac count (similar t o Friend & Abbotliri986L 
IPauldrach et alj[l986l and lKudritzki et al Jl I989f) . During these two 
first steps of the series, the ionisation and excitation calculation is 
further simplified by setting all recombination fractions C,i = 1 and 
by assuming a Boltzmann excitation formula [i.e. W = 1 in Equa¬ 
tion GoJj. Finally, in the last stage of the series, the assumption 
of an unattenuated radiation field will be dropped. In this calcula¬ 
tion, we then use the full ionisation and excitation description as 
outlined in Section [X71 including recombination fractions and the 
geometric dilution factor in the excitation balance. This leads to a 
mild variation of the degree of ionisation in the wind. When ap¬ 
plicable, the results of the Mcrh calculations will be confronted 
with the predictions of the CAK and Mcak theory. These compar¬ 
isons are performed on the basis of the stationary wind structure 
found in all calculations. With this in mind, we further reduce the 
computational complexity in all Mcrh simulations by following 
all Monte Carlo packets until they escape the wind during each 
time step. This is equivalent to assuming a light propagation speed 
much larger than the fluid flow and constitutes a common strat¬ 
egy in Mo nte Carlo radiative tra nsfer approache s, for example in 
PYTHON dLong & Kniggel r20021 or TARDIS iKerzendorf & Sinj 
120141) . Adopting this technique absolves us of the need to invest 
computational power into the initial build-up of the radiation field. 
As long as we are solely interested in the final stationary wind state, 
this simplification is justified. For future calculations, which aim at 
investigating the dynamical behaviour of the wind structure, this 
modification of the Monte Carlo radiative transfer scheme will be 
dropped. 


4.4 General parameters 

In the Mcrh simulations, the evolution of the wind material is 
followed until a steady-state structure emerges. All calculations 
are carried out on a non-uniform spherical mesh with 512 cells, 
which span the region between 1 77* and 10 77*. The cell size in¬ 
creases outwards from 1.76 x 10~ 3 77* at the stellar surface to 
1.73 x 10 _1 77* at the outer edge of the domain. During each time 
step, the incident radiation field is discretised by 5 x 10 4 packets, 
which sample the black-body spectrum in the wavelength range be¬ 
tween A m in = 228 A and A max = 22800 A. The lower edge corre¬ 
sponds to the ionisation edge of He II. Any radiation more energetic 
than this threshold is assumed to be rapidly removed by bound-free 
absorptions by helium and consequently not included in our consid¬ 
eration (see|Simj|2004j, for a similar strategy). For the initial wind 
state in the MCRH calculations, we adopt a structure that is simi¬ 
lar to the Cak/Mcak solution but scaled up or down. By experi¬ 
menting with other initial wind configurations, we have explicitly 
verified that the final stationary wind solution, which is found in 
the calculations, is insensitive to the details of the initial state. Fur¬ 
thermore, we have also ensured that enough Monte Carlo packets 
are used in the simulations. Increasing the number does not change 
the overall shape of the radiative acceleration but only reduces the 
strength of the stochastic fluctuations (see Figure[4jl. 


4.6 Mcrh Cak/Mcak module 

To better judge the outcome of the Mcrh-Cak/Mcak compar¬ 
ison, we include an additional set of numerical calculations per¬ 
formed with our approach. Instead of relying on the Monte Carlo 
radiative transfer procedure to determine the line-driving force, we 
incorporate the analytic power-law expressions of the Cak/Mcak 
approach. In each simulation cycle, after the fluid dynamical and 
central gravity substeps, the instantaneous value of the dimension¬ 
less optical depth parameter t is determined in each cell and the 
force multiplier M (t) calculated accordingly. In this procedure, the 
central star may be treated either as a point source [see Equation 
{Al l, or its finite extent may be taken into account [see Equations 
< |A8] ) and jA9l l. To determine the required instantaneous velocity 
gradient in each cell, we rely on the algorithm bv lFornbergl ( t 198a) 
to construct finite differences on arbitrarily spaced grids. We will 
refer to all numerical radiation hydrodynamical calculations per¬ 
formed with this module as Cak-Rh and Mcak-Rh respectively. 

4.7 Results 

4.7.1 Calculation in the point-source limit 

In the first set of Mcrh calculations, the central star is approxi¬ 
mated by a point source and the radiation field is assumed to be 
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unattenuated. Within the Monte Carlo approach, both simplifica¬ 
tions are realised by injecting the packets only on radial trajecto¬ 
ries, i.e. replacing the random number experiment © by n = 1, 
and by not changing the propagation direction in interaction events, 
i.e. i-iQ fiQ. 

In Figure [3] the final stationary wind state, which establishes 
in the Mcrh calculation, is shown in terms of the density, velocity 
and mass-loss rate. The comparison with the results of the Cak- 
Rh simulation and the analytic predictions according to the Cak 
theory shows a very good agreement. This positive finding is a first 
indication for the accuracy and utility of our Monte Carlo-based 
scheme to address the line-driving problem in hot-star winds self- 
consistently. 


4.7.2 Including the finite-cone effect 

During the second stage of the wind simulation series, the point- 
source approximation is dropped and the finite extent of the central 
star is taken into acco unt, analogously to the Mc ak approach of 
IPauldrach et alj dl986h and IKudritzki et all dl989h . In the Mcrh 
calculations during this stage, we allow Monte Carlo packets to 
propagate on non-radial rays as well by sampling the initial direc¬ 
tion according to Equation ©. 

As in the point-source calculations, the stationary wind struc¬ 
ture obtained with MCRH is compared with the analytic predic¬ 
tions, now according to the Mcak theory as outlined in Appendix 
lAl Again, we include the numerical results calculated with the 
Mcak-Rh version of our scheme. When accounting for non-radial 
photon propagation paths, the Monte Carlo-based results agree 
again very well with the analytic predictions and the Mcak-Rh 
calculation, as illustrated by Figure [3] 

As ex pected from numerous previous studies , mos t no¬ 
tably from iFriend & Abbott! dl986h , IPauldrach et alj dl986h and 
iKudritzki et aid d 19891) , the inclusion of the finite-cone effect leads 
to higher wind velocities, in our case by a factor of 2.5 in terms of 
the velocity at r = 10 f?* compared to the point-source case. At 
the same time, the mass-loss rate drops by a factor of 0.5. Figure 
[4]illustrates the difference in the radial dependence of the radiative 
acceleration, which is responsible for the change in the structure of 
the wind flow. 


4.7.3 Full inclusion of scatterings 

After having established the basic applicability of our approach to 
the line-driving problem in the first stages of the simulation series, 
we take another step towards a realistic description of the radiation 
field in hot-star winds by dropping the unattenuation approxima¬ 
tion and accounting for the full scattering process in interactions. 
Now the entire procedure described in Section [3~4l is performed: in 
particular, the emergent propagation directions in atomic line inter¬ 
actions are drawn according to Equation GS- Accounting fully for 
the multiple scattering phenomenon, these calculations reach be¬ 
yond the capabilities of the basic CAK and Mcak approaches (as 
outlined in Appendix[A]l. We emphasise, that the key consequence 
of multiple scattering lies in the capability of line interactions to 
lengthen the photon propagation trajectory. By this process, pho¬ 
tons may potentially exert a stronger acceleration onto the wind ma¬ 
terial, compared to the unattenuated case in which they may also in¬ 
teract multiple times but their trajectories are always straight lines. 
Once the unattenuation of the radiation is dropped in the Mcrh 
simulations, some packets may potentially be back-scattered onto 



Figure 4. The radial dependence of the line-driving force in the MCRH 
calculations after a stationary state has established. The colour-coding is 
the same as in Figure Q] For the results shown by the weak solid lines, 
5 X 10 4 packets where used to discretise the radiation field emitted by the 
central star in each simulation cycle. For comparison, the full solid lines 
show the corresponding acceleration, when the 10 6 packets are used once a 
stationary state has emerged. Notice the clear decrease in the Monte Carlo 
noise but the same radial dependence of the radiative acceleration as in the 
case with fewer packets. 


the stellar disc. To counteract the luminosity loss induced by this 
process, we follow iLucv & Abbott! <[ 1993 T and scale the packet en¬ 
ergies by a constant factor, which ensures that in each time step, net 
radiative energy amounting to the luminosity L* is streaming into 
the wind. This process may be interpreted as a colour correction of 
the stellar spectral energy distribution (see lLuc vl ll999bl for a sim¬ 
ilar strategy in the context of calculating synthetic observables for 
supernovae). 

Accounting both for the finite-cone effect and the full scat¬ 
tering procedure, we again determine the stationary wind structure 
with Mcrh. The resulting wind velocity, the density stratification, 
and the mass-loss rate are included in the summary plot of Figure 
[3] Compared to the unattenuated calculations with the finite-cone 
effect, we find a slightly slower wind and, as shown in Figure [4] 
a weaker radiative acceleration. At first glance, this seems to con¬ 
tradict the statement about the effect of multiple scattering in the 
introductory part of this section. But one has to bear in mind that in 
the Cak/Mcak description of the line-driving problem, in each in¬ 
teraction the full photon/packet momentum is transferred onto the 
wind material. This is comparable to a purely absorbing medium, 
with the important difference that the photon trajectory is not termi¬ 
nated in the Cak/Mcak description. Instead, the same photon may 
still contribute to the acceleration in regions of the wind located at 
larger radii. In the full scattering case, in which the re-emission of 
the line-interaction is taken into account, no momentum is trans¬ 
ferred onto the wind on a Cak/Mcak like trajectory since it in¬ 
volves forward-scatterings only. Thus, the momentum transfer in 
the Cak/Mcak case may generally be overestimated and only the 
photons that are on trajectories that have been significantly length¬ 
ened due to many successive scatterings may contribute compara¬ 
bly to the Cak/Mcak description. This phenomenon is illustrated 
in Figure [5] A general reduction of the line-driving force once the 
unattenuation ass umption is dropped has already been found in pre¬ 
vious studies, e.e. lAbbott & Lucvl lil985l ). 
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MCRH results 

Cak-Rh/Mcak-Rh results 


Cak/Mcak predictions 
point source, unattenuation 


finite cone, unattenuation 
finite cone, full scattering 
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Figure 3. Results of the stellar wind simulation series in terms of the final stationary wind velocity (left panel), density (central panel) and mass-loss rate (right 
panel). The colour-coding reflects the different stages of the simulation series. The calculations which are based on the point-source approximation and use 
the unattenuation of the radiation field are shown in blue. Red lines corresponds to calculations, in which the finite-cone effect is included and green to those 
which also include the full scattering procedure. All MCRH results are presented by solid lines. Where applicable, the analytic predictions according to the 
Cak and the MCAK theory are included as dotted lines. In these cases, also the results obtained with Cak-Rh/Mcak-Rh are given as an additional reference 
(dashed lines). 
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Figure 6. Detailed summary of the contribution of the different elements to the line-driving force in the full-scattering simulation. In the left panel, the 
composition of the line-driving force averaged over the entire simulation box is displayed. Here, also the contributions of the different ionisation stages are 
shown. The right panel shows how this composition varies throughout the wind. For this illustration, the contributions to the line-driving for of all ions of one 
element are combined and then colour-coded. Moreover, a binning with a width of ten grid cells has been used. 


4.7.4 Additional diagnostics 

By recording the details of all interaction events performed by the 
Monte Carlo packets in the final calculation of the simulation se¬ 
ries, the origin of the radiative acceleration can be studied. The 
contributions of the different elements and ionisation stages, aver¬ 
aged over the entire computational domain, are illustrated in Figure 
[6] This highlights that the line-driving force mainly derives from 
interactions with lines of iron, nickel some intermediate mass ele¬ 
ments and the CNO group. The relative importance of the differ¬ 
ent contributions, however, changes throughout the wind, as shown 
in the right panel of the same figure. In our simulation, lines of 
iron group elements mostly contribute in the inner part of the wind, 
close to the photosphere. Further out, the intermediate mass ele¬ 
ments grow in importance. This finding is compatible with the in¬ 
vestigation of Ivink et alj dl999h . in which the importance of iron 
for the radiative acceleration in the lower parts of the wind has been 


highlighted in the context of the bi-stability jump. We stress, how¬ 
ever, that due to the simplified treatment of ionisation and excita¬ 
tion in our simulations, the results presented in Figure[6]should not 
be over-interpreted, but viewed as an illustration of the diagnostic 
possibilities of our approach. 

By recording the interaction histories of all packets we can 
also investigate the importance of multiple scattering in our simu¬ 
lation. FigureQshows the number of interactions performed by the 
packets that escaped through the outer edge of the computational 
domain. In our particular simulation, most packets never explic¬ 
itly interacted on their way out. They may have still contributed to 
the line-driving force as long as they came into resonance with at 
least one line transition (see discussion in Section [T6l >. The major¬ 
ity of those that interacted, however, performed multiple scattering 
events. 
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Figure 5. Illustration of the total radiation momentum transfer rate exerted 
by the Monte Carlo packets along their propagation paths in terms of a two- 
dimensional histogram. The length of the trajectories, l, is normalised to 
the radial extent of the domain, l r . The upper two panels show the accumu¬ 
lated momentum in the point-source (left) and the finite-cone (right) MCRH 
calculation respectively. The lower panels correspond to the MCRH wind 
simulation with the full scattering procedure. The left panel only includes 
packets that are ultimately backscattered onto the central star. Thus, trajec¬ 
tories shorter than l r are encountered. In the lower right panel only packets 
which ultimately escape through the outer boundary of the computational 
domain are shown. 
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Figure 7. Histogram of the number of interactions which Monte Carlo pack¬ 
ets that ultimately escape through the outer edge of the simulation box per¬ 
form. As displayed in the inset, the majority of the escaping packets never 
experience an interaction. However, most of the interacting packets perform 
multiple scatterings. 
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5 COMPARISON WITH DETAILED CALCULATIONS 

Having established the basic performance of our method under sim¬ 
plified conditions, we now compare our approach with modern, 
sophisticated techniques for determining the struc ture of hot-star 
winds , in particular with the approach developed bv lMuller & Vinkl 
(2008|, MV08 hereafter). This method constitutes an advance¬ 
ment o f the original approach of Ivink et al.l i ll 9991) and I Vink et alj 
(2000). A Monte Carlo radiative transfer calculation is used to de¬ 
termine the local line-driving force for a given static wind structure. 
Here a velocity structure that is more general and flexible than the 


Table 4. Simulation parameters for the calculation of the wind from the 
05-V main sequence star Ic.f. lMtiller & Vinkl2008l) . 


Parameter 

Value 

M* 

40 Mq 

log T+/L q 

5.5 

TeS 

4 x 10 4 K 

Ne/W 

10 14 cm” 3 

P0 

10 -12 gem -3 


/3-type law is used. According to the reconstructed line accelera¬ 
tion, the mass-loss rate and the parameters of the wind velocity law 
are iteratively updated until a converged wind structure has been 
found. 


5.1 Parameters 

MV08 test their approach by predicting the wind structure of an 
05-V main sequence star. For direct comparison, we adopt param¬ 
eters in our simulation to match their calculation. In Table [4] all 
quantities which have been changed with respect to the calculations 
presented in Section[4]are listed. These choices imply an Eddington 
factor of T e = 0.210, which is very close to the value quoted by 
MV08. 

The comparison with the Mu ller & Vinkl J2008l> technique will 
be based on MCRH calculations that incorporate all techniques de¬ 
scribed in Section]!] Since no Cak or MCAK-like simulations are 
performed in this context, the CAK fitting procedure of Section H31 
does not have to be repeated for the current stellar parameters. 


5.2 Results 

With the stellar parameters listed in Table[4]a full scattering MCRH 
simulation, similar to the calculation presented in Section 14.7.31 
is performed to determine the structure of the wind of the 05-V 
main sequence star. The result is shown in Figure [8] in terms of the 
stationary velocity and mass-loss rate and compared to the structure 
found by MV08. To obtain the comparison MV08 wind velocity, 
we use the approximate expression, Equation 39 in MV08, which 
is only strictly valid in the supersonic wind regime. 

As seen in Figure [8] both approaches predict winds that are 
quite similar in shape. However, the velocity structure found by 
MCRH accelerates quite quickly towards the terminal wind speed, 
whereas the wind velocity law predicted by MV08 approaches its 
final value in a gentler fashion. The two winds deviate slightly 
in absolute values: the MCRH wind reaches a velocity of u = 
3065kms _1 at r = 10 f?*, while the MV08 solution lies at u = 
2719 km s -1 . The mass-loss rates of the winds agree on a similar, 
namely a ten percent level, with M = 7.97 x 10 -7 Mq yr _1 being 
obtained in the MCRH calculation and M = 8.99 x 10 -7 Mq yr _1 
found by MV08. 

These minor discrepancies may partly be related to the dif¬ 
ferent solution strategies followed in the two approaches. In our 
method, a radiation hydrodynamical calculation is relaxed to a 
steady state solution without any major restrictions on the shape 
of the line acceleration and velocity structure. The MV08 approach 
relies on a parameterisation of the line driving force and in turn 
of the velocity law, thus allowing only wind structures of a certain 
family. Differences between the two calculations will also arise be¬ 
cause of the very simplistic treatment of ionisation and excitation in 
our method. Overall, however, the level of agreement between the 
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Figure 8. Stationary wind structure for the 05-V main sequence star deter¬ 
mined with Mcrh (green). The predictions according to the MV08 tech¬ 
nique are shown as a comparison (grey). 


packets) were executed per cycle. By comparison, only 0.015 s are 
required and 1.3MFlops executed per cycle if the Cak-Rh ver¬ 
sion is used instead. We emphasise, however, that all calculations 
presented in this work were carried out serially on a single pro¬ 
cessing core. We stress, that the performance of the Monte Carlo 
scheme may be significantly improved by devising a parallelisa¬ 
tion scheme and distributing the workload onto many cores. For all 
calculations, the executable was produced with the gcc compiler, 
version 4.7.3, using the optimisation level -03. 


6.2 Influence of Monte Carlo noise 


The probabilistic nature of Monte Carlo techniques leads to the in- 
evitable introduction of s tochastic fluctuations (see, for example, 
ICarter & Cashwelj j 1975) for a disc ussion of Monte Carlo errors). 
However, Noebauer et al. H20li ) and lRoth & Kasenll2015l) demon- 
strated that, with appropriate reconstruction techniques, the fluctu¬ 
ations can be controlled and Monte Carlo-based radiation hydrody- 
namical simulations can indeed be performed. The simulations per¬ 
formed in this work and presented in the previous section illustrate 
that the Monte Carlo noise, which is clearly visible in the recon¬ 
structed radiative acceleration (c.f. Figure 0, does not prohibit the 
use of Monte Carlo methods to study the line-driving problem. In 
fact, the obtained velocity and density structure (shown in Figure[3]l 
is remarkably smooth. Due to their pure stochastic nature the fluc¬ 
tuations cancel in an average sense (but see also remarks in Section 

IQt . 


6.3 Full non-LTE treatment 


two calculations is good and lends credence to their complementary 
use in the study of stellar winds. 


6 DISCUSSION 

The results of the wind calculations and the diagnostic capabilities 
presented in the previous sections, in particular the good agreement 
between the MCRH results and the Cak/Mcak predictions and the 
comparison to the MV08 technique, clearly demonstrate the utility 
and accuracy of our Monte Carlo-based radiation hydrodynamical 
approach for solving the line-driving problem in hot-star winds. In 
the following, we briefly discuss some particular features to our 
approach and also comment on some of its current limitations. We 
also sketch possible alterations of our approach aimed at alleviating 
these shortcomings and discuss potential future applications. 


6.1 Numerical performance 

In general, Monte Carlo radiative transfer techniques are rather 
computationally demanding, since many Monte Carlo quanta have 
to be processed to reach the desired statistical fidelity. This draw¬ 
back is typically balanced by the very favourable parallelisation 
properties of such calculations (see discussion in section IQl . In 
our case, the Monte Carlo-based radiation hydrodynamical simu¬ 
lations are indeed costly. In the full scattering simulation, the fi¬ 
nal stationary state established after about 4800 simulation cycles. 
On of these takes roughly 14 s on a Intel Xeon E5520 processor 
if 5 x 10 4 packets describe the incident radiation field. This in¬ 
creases to about 210 s once 10 6 packets are used. In these calcula¬ 
tions about 10 GFlops (for 5 x 10 4 packets) and 160 GFlops (10 6 


The calculations presented here assumed either LTE or relied on 
the approximate non-LTE prescription presented in Section 13.71 
Non-LTE effects, however , play an important role in stellar winds 
(see, e.g., IPauldracjl 19871 : IPauldrach et alj[l994l : IPuls et alll2005l . 
for non-LTE calculations of hot-star winds) and should be included 
in future calculations. 

In the general non-LTE situation, the level popula t ions o f ex¬ 
cited states are influenced by the radiation field. iLuc vl d2002l) and 
Luc 3 Hqo 1) describes a simple and elegant procedure to recon¬ 
struct these radiative rates from a Monte Carlo simulation. In a first 
step of refining the physical detail in the interaction treatment of our 
approach, this procedure could be adopted to determine a non-LTE 
excitation and ionisation balance during the Monte Carlo radiative 
transfer step. This strategy has already been followed in p revious 
Monte Carlo-base d studies, for example b y lsim et al] d2005h . using 
the PYTHON code dLong & K nigge 20()2|) an d in non-LTE radiative 
transfe r calculations in Keplerian discs by ICarc iofi & Bio rkmanl 

d 2006 t) . 


To go beyond the pure resonant scattering approximation in 
th e line- i nteract ion tr eatment, either the simple branching scheme 
of lLucvl dl999bl ) (see lMazzaiill200(]| . for Monte Carlo-based radia¬ 
tive transfer calculations in Type la ejecta using this scheme) may 
be incorporated or the so-called Macro-Atom formalism introduced 
by L uc\d (2 002j , 120031 ) may be used. The latter fully accounts for 
all processes that may excite or de-excite line-transitions, includ¬ 
ing collisional processes. We emphasise, that the utility and ac¬ 
curacy of thes e techniques haf^been carefully examined and es¬ 
tablished dLucy||2003 . 120031 , 120051) . Moreover, they have already 
been incorporated into a number of Monte Carlo radia tive trans¬ 
fer frameworks, including Artis i Kromer_&Mm||200S ), PYTHON 
dSim et al.ll2005t) and TARDIS dKerzendorf & Sirnl2014 ). 
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6.4 Multidimensionality and non-Sobolev treatment 


Observations and also theoretical investigations indicate that line- 
driven winds are neither stationary nor smooth, but feature avari- 
able and inhomogeneous structure (see, e.g., overview bv lPuls et alj 
120081 ). These findings advocate a dynamical and fully multidimen¬ 
sional treatment of the problem. Monte Carlo techniques are suited 
for this task since the formalism readily generalises (conceptually) 
to arbitrary geometrical configurations and also exhibits excellent 
scaling properties on parallel processors - a desirable feature when 
facing the increased numerical workload of multidimensional cal¬ 
culations (see, for example, comments by Kasen et al. 2 0061 and 
|BaesetaI] 20H, as well as the scaling tests bv iRoth & Kasenl 20151 
and Harries! 2015t) . With this in mind an obvious extension of our 
method comprises a generalisation to multidimensional geometries 
to address inhomogeneous line-driven outflows. In this context, 
some work should be invested into devising an efficient interfac¬ 
ing of the standard parallelisation strategies of Monte Carlo and 
fluid dynamical techniques (see e.g. discussion in lBaes et alj|201l1 


and developments by lHarriesll20 1 51 . ). 

When investigating multidimensional line-driven mass out¬ 
flows, a generalisation of the line interaction treatment that goes 
beyond the Sobolev approximation used here should also be con- 
sidered: this would be necessary to stud y the line-driving instability 
dLucv & Solomonl 1970i: lo woe kil 1 9941) which should occur in line- 
driven outflows and may play a part in understan din g the clump - 
ing mechanism (see discussion in lPuls et alJI 2008! and lvinkll2015h . 


Non-Sobolev Monte Carlo radiative transf er sc hemes have already 
been developed and u sed, for example bv lKnigge et al J i ll9951) and 
iKusterer et alj 120141) in the context of winds of cataclysmic vari¬ 
ables. However, the computational costs of such schemes are sig¬ 
nificantly higher than their Sobolev counterparts. Moreover, the in¬ 
herent Monte Carlo noise could potentially be more problematic 
than described in 16.21 since small perturbations in the line-driven 
wind may, in principle, self amplify in non-Sobolev calculations. 
Further investigation is required to asses the relevance of this po¬ 
tential caveat. 


Mcak theory. We demonstrated that our method can also go be¬ 
yond the capabilities of Cak and Mcak by dealing with an atten¬ 
uated radiation field and the effects of multiple scattering. For the 
particular physical conditions investigated in Section [4] these ef¬ 
fects lead to a reduction of the line-driving force and thus a slower 
wind velocity structure compared to the Mcak predictions. Finally, 
we compared results obtained with our approach to those computed 
by MV08 for an 05-V main sequence star and found good agree¬ 
ment, given the difference in approach and simplifications made 
here. 


The successful outcome of our wind simulations demonstrates 
the possibility to use a Monte Carlo-based radiation hydrody- 
namical approach to model line-driven mass outflows. Thus, this 
approach holds promise for detailed multidimensional and self- 
consistent studies of inhomogeneous winds, including multi-line 
effects. With a fully multidimensional version of our approach, 
line-driven outflows in systems other than hot-star winds may 
be addressed as well. In particular, the winds emanating from 
accretion discs in ca taclysmic variables (e.g. IProga et aH 1 


NoebaueretaL Il20ld) or active galactic nuclei (e.g. IProga et alj 
2000l : IProga & Kallmanll2004l : iHigginbottom et all |20 1 3l> may be 


investigated self-consistently and in great detail with such a Monte 
Carlo-based approach. 
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7 CONCLUSIONS 

In this work we have introduced a new approach to solve line- 
driven stellar winds self-consistently by using a Monte Carlo-based 
radiation hydrodynamical approach. The key feature of this tech¬ 
nique lies in the reliance on a Monte Carlo radiative scheme. This 
technique offers a number of advantages when dealing with com¬ 
plex interaction physics, in particular when multi-line effects play 
a role. Moreover, this Monte Carlo-based technique readily gener¬ 
alises to multidimensional geometries, a very advantageous feature 
for potential studies of inhomogeneous outflows. 

Establishing the utility and accuracy of the introduced ap¬ 
proach in solving the line-driving problem was the main focus of 
this work. Consequently, we designed the calculations to capture 
the essence of the line-driving problem, but adopted a number of 
simplifications. These do not interfere with the general applicabil¬ 
ity of our method to solve the local line-driving problem, but re¬ 
duce the computational complexity. Most importantly, we adopted 
the Sobolev approximation and a simple and approximate non-LTE 
treatment to determine the excitation and ionisation balance. 

Using our scheme, we successfully solved for the stationary 
structure of one-dimensional, spherically symmetric hot-star winds 
achieving good agreement with the predictions of the CAK and 


APPENDIX A: CAK AND MCAK APPROAC HE S 

In the following, the relevant expressions of the CAK theory are 
listed in the form in which they are used in the current study. In 
general, the line-driving force is expressed as a multiple of the ac¬ 
celeration due to electron scattering 

ref r 

flline = M (t ) ^2,, ’ (A1 ^ 

assuming a reference specific interaction cross section al el . 
Throughout this work, <rf f = 0.3 cm 2 g -1 is used. The force mul¬ 
tiplier depends on the dimensionless optical-depth t 

t = a T e el pu th > (A2) 

which involv es the thermal velocity of the wind material (c.f. 
lAbbodll982h 

/ 2 k B T eS 

u±h = \ -, (A3) 

V m H 

and is approximated by the power-law 

M(i) = kt~ a (A4) 
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if the cent ral st ar is assumed to radiate as a point-source (again fol- 
lowing lcastor et aDl 19751 ). It should be noted that, throughout this 
work, we neglect any influence of changes in the ionisation, which 
is equivalent to setting S = 0 in the force multiplier formulation 
proposed bv lAbbottl <t 1982h . 

In the point-source limit, the solution to the momentum equa¬ 
tion CD, with the Cak line-driving force, results in the wind veloc¬ 
ity law 


u(r) = Uoot/l - 


R* 


(A5) 


with the terminal speed being a multiple of the local escape speed 
from the photosphere, u eS c: 


U oo — 


a 


~ Uesc • 


1 — a 

The constant mass-loss rate of the wind is given by 

Mcak = — ( — GM.il - T< X 

^th \ ® e ref 


(A 6 ) 


k -u 

c 


x a( 1 — a)" 


(A7) 


Once the finite extent of the star is taken into account, the force 
multiplier is modified by a finite-cone correction factor 


Mpc(t) = Dpc{t)McAK{t), (A 8 ) 


Dpc(t) 

which involves 


(l + a ) a+1 - (l + apl ) a+1 
(1 — pl)(a + 1 )<t (1 + er)“ 


cr 


■ drt 
u dr 


M* = 



(A9) 


(A10) 

(All) 


We follow lKudritzki et al J j 1989h and predict the wind structure ac¬ 
cording to their approximate analytic solution technique. In partic¬ 
ular, we adopt their proposed approach for the case of a “frozen-in” 
ionisation state (5 = 0, c.f. Ku dritzki et all 19891 section 4.1). The 
wind velocity is assumed to be very close to a /3-type law and the 
mass-loss rate decreases with respect to the original CAK case ac¬ 
cording to 


M = 


1 + a 


Mcak ■ 


(A12) 


The wind velocity in this approach follows from performing the 
integration 


In these expressions, the inverse of the radial distance relative to 
the photosphere, v = R+/r, is used. Throu ghou t this work, the 
Mcak equations are solved for /3 = 0.8 (see lPauldrach et ali 19861 : 
Kudritzki et alD l 19891 for a motivation of this value). 


APPENDIX B: RADIATION HYDRODYNAMICS 
EQUATIONS 


We briefly review the radiation hydrod ynamical equations on w hich 
our numerical scheme is based (c.f. Miha las & Mihalasll 19841) . In 
general, these equations describe the conservation of mass, momen¬ 
tum and energy. Since we assume in this work, that the wind out¬ 
flow remains isothermal at all times, we only consider the first two 
of these equations and use the isothermal equation of state to re¬ 
late fluid density and thermodynamic pressure. Cast into a pseudo- 
Lagrangian form by using the substantial derivative 


D_ _ d_ _d_ 

Df df + dr ’ 


(Bl) 


this radiation hydrodynamical pro blem in one-dimensional spheri- 
cal symmetry is described by (c.f. lMihalas & MihalaslfT984h 


D 1 d . 2 n 
P + P— ~ ( r w ) = °. 


Df 


r 2 dr' 

D d j 

PTT+ U + 1Z P = P9 + G , 


Df 


dr 


T~) 2 

P = d;« 


(B2) 

(B3) 

(B4) 


The fluid density p and its velocity u appear together with the 
thermodynamic pressure P, the isothermal sound speed di ao and 
a static external gravitational field g. The transfer of momentum 
between the fluid and the radiation field is captured by the radia¬ 
tion force, which ma y be determined from the first moment of the 
transfer equation (c.f. lMihalas & MihalaJl984i) : 

c\_ roo n 1 

G 1 = — j d V J d//(x / - v)p- (B 5 ) 

Here, the description of the radiation field by the specific intensity 
I is used and the material functions opacity x and emissivity 77 , 
describing the absorption and emission of radiative energy, appear. 
The integration is performed with respect to the entire frequency 
spectrum and to all possible values for the cosine of the propagation 
direction, p. 
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